---
title: "Nighttime temperature and human sleep loss in a changing climate"
subtitle: "Replication code"
author: Nick Obradovich^[Corresponding author.] ^[Belfer Center for Science and International Affairs, Kennedy School of Government, Harvard University.] ^[Media Lab, Massachusetts Institute of Technology.], Robyn Migliorini^[ Department of Psychiatry, Massachusetts General Hospital] ^[Harvard Medical School] ^[San Diego State University/University of California San Diego Joint Doctoral Program in Clinical Psychology.], Sara C. Mednick^[Department of Psychology, University of California Riverside.], James H. Fowler^[Departments of Political Science and Medicine, University of California San Diego.]
date: ""
output: 
  html_document:
    fig_caption: true
---

```{r load, echo=T, message=F, warning=F, include=T}
knitr::opts_chunk$set(fig.path='figures/', warning=F, message=F, fig.retina=1, cache=T, autodep=T, echo=F)

rm(list = ls())

library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
#library(kfigr) #figure referencing for markdown
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots

# data and functions
load(file='./climate_sleep.Rdata')

# set cores for parallel processing
registerDoMC(20)

# run from command line
# Rscript -e "rmarkdown::render('./main_text_replication_code.Rmd')"
```

```{r figure1, echo=TRUE, fig.width=10, fig.height=6.5, dev='jpeg', dpi=1200, fig.cap="**Main text Figure 1**"}
# main model
model <- felm(qlrest2 ~  avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
              | time + mmsa:season
              | 0 
              | mmsa + time, data=brfss, exactDOF=TRUE, psdef=FALSE)
summary(model)

# qlrest2 graphs

# tmin
brfss <- brfss[, tmin.dev.round := round_any(avg.tmin.dev, 0.5)] # round to bins of 0.5 degree
tmin.binned <- brfss[, .(avg.qlrest2 = mean(qlrest2, na.rm=T),
                         high.qlrest2 = mean(qlrest2, na.rm=T) + sem(qlrest2),
                         low.qlrest2 = mean(qlrest2, na.rm=T) - sem(qlrest2)), 
                     by=.(tmin.dev.round)] # average and sem dv over each bin

tmin <- ggplot(data=tmin.binned[abs(tmin.dev.round) <= 5,], aes(x=tmin.dev.round, y=avg.qlrest2)) + 
            geom_smooth(data=brfss[abs(avg.tmin.dev) <= 5,], aes(x=avg.tmin.dev, y=qlrest2), fill="goldenrod1", color="firebrick", size=1.75, alpha=1) + 
            geom_point(color="firebrick", size=3.5) + 
            ggtitle('a') +
            ylab("Mean Monthly Nights of Insufficient Sleep") +
            xlab("Mean Monthly Nighttime Temp.\n(Anomaly in \u2103)") +
            coord_cartesian(xlim=c(-5.2, 5.2), ylim=c(6.95, 8.35), expand=FALSE) +
            theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
                plot.title = element_text(hjust=0.025),
                axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
                axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
                axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
                axis.text.x = element_text(size=18, face="bold"),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                panel.background = element_blank(),
                axis.line.x = element_line(),
                axis.line.y = element_line(),
                legend.title=element_blank()
              )

# tmin histogram
tmin.hist <- ggplot(data = brfss, aes(x = avg.tmin.dev)) +
                  geom_histogram(aes(y=..density..), fill='firebrick', color='firebrick4', alpha=1) +
                  ylab("Density") +
                  xlab("Mean Monthly Nighttime Temp.\n(Anomaly in \u2103)") +
                  coord_cartesian(xlim = c(-10, 10), ylim = c(-0.0025,.31), expand=FALSE) +
                  scale_x_continuous(breaks=c(-5,-2.5,0,2.5,5)) +
                  ggtitle("b") +
                  theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
                        plot.title = element_text(hjust=0.025),
                        axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
                        axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
                        axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
                        axis.text.x = element_text(size=18, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        legend.title=element_blank()
                  )

grid.arrange(tmin, tmin.hist, ncol = 2)
```

```{r figure2, echo=TRUE, fig.width=16, fig.height=6.5, dev='jpeg', dpi=1200, fig.cap="**Main text Figure 2**"}
# qlrest2 by season, income, and age coefficient graphs

# seasons
spring <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
               | mmsa + time
               | 0 
               | mmsa + time
               , data=brfss[season=="spring",], exactDOF=TRUE, psdef=FALSE)
summary(spring)

summer <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud  
               | mmsa + time 
               | 0 
               | mmsa + time
               , data=brfss[season=="summer",], exactDOF=TRUE, psdef=FALSE)
summary(summer)

fall <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud  
             | mmsa + time
             | 0 
             | mmsa + time
             , data=brfss[season=="fall",], exactDOF=TRUE, psdef=FALSE)
summary(fall)

winter <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud  
               | mmsa + time
               | 0 
               | mmsa + time
               , data=brfss[season=="winter",], exactDOF=TRUE, psdef=FALSE)
summary(winter)

# get coefficients of variables
beta.spring = spring$coefficients[row.names(spring$coefficients)=="avg.tmin.dev"]
beta.summer = summer$coefficients[row.names(summer$coefficients)=="avg.tmin.dev"]
beta.fall = fall$coefficients[row.names(fall$coefficients)=="avg.tmin.dev"]
beta.winter = winter$coefficients[row.names(winter$coefficients)=="avg.tmin.dev"]

# cluster robust standard errors
se.spring = spring$cse[names(spring$cse)=="avg.tmin.dev"]
se.summer = summer$cse[names(summer$cse)=="avg.tmin.dev"]
se.fall = fall$cse[names(fall$cse)=="avg.tmin.dev"]
se.winter = winter$cse[names(winter$cse)=="avg.tmin.dev"]

# upper and lower confidence bounds
z.score.spring = qnorm(1 - ((1 - 0.68)/2))
upper.spring = beta.spring + z.score.spring*se.spring
lower.spring = beta.spring - z.score.spring*se.spring
z.score.summer = qnorm(1 - ((1 - 0.68)/2))
upper.summer = beta.summer + z.score.summer*se.summer
lower.summer = beta.summer - z.score.summer*se.summer
z.score.fall = qnorm(1 - ((1 - 0.68)/2))
upper.fall = beta.fall + z.score.fall*se.fall
lower.fall = beta.fall - z.score.fall*se.fall
z.score.winter = qnorm(1 - ((1 - 0.68)/2))
upper.winter = beta.winter + z.score.winter*se.winter
lower.winter = beta.winter - z.score.winter*se.winter

d <- as.data.table(cbind(rbind(beta.spring,beta.summer,beta.fall,beta.winter),rbind(upper.spring,upper.summer,upper.fall,upper.winter),rbind(lower.spring,lower.summer,lower.fall,lower.winter),rbind("Spring","Summer","Fall","Winter")))
setnames(d, c('beta', 'upper_bound', 'lower_bound', 'variables'))
d <- d[, beta := as.numeric(beta)*100] # multiply by 100 to make more interpretable (effect per 100 individuals)
d <- d[, upper_bound := as.numeric(upper_bound)*100]
d <- d[, lower_bound := as.numeric(lower_bound)*100]
d <- d[, variables := factor(variables, levels=c("Spring","Summer","Fall","Winter"))]

season <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        ggtitle("a") +
        geom_bar(stat = "identity", aes(width = 0.80), color='transparent') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("gray70","firebrick","gray70","gray70"),breaks=c("Spring","Summer","Fall","Winter")) +
        coord_cartesian(ylim = c(-1,10), expand=TRUE) +
        ylab("Marginal Effect of \u2103 Anomaly\non Monthly Nights of Insufficient Sleep\n(per 100 individuals)") +
        xlab("Season") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=16,angle=0, vjust=.5, face="bold") ,
              axis.text.x = element_text(size=15, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.y = element_line(),
              axis.line.x = element_line(),
              legend.position="none"
        )   

# income
low.income <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
                   | time + mmsa:season 
                   | 0 
                   | mmsa + time
                   , data=brfss[income2 %in% 1:6,], exactDOF=TRUE, psdef=FALSE)
summary(low.income)

high.income <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
                    | time + mmsa:season
                    | 0 
                    | mmsa + time
                    , data=brfss[!(income2 %in% 1:6) & !is.na(income2),], exactDOF=TRUE, psdef=FALSE)
summary(high.income)

# get coefficients of variables
beta.low = low.income$coefficients[row.names(low.income$coefficients)=="avg.tmin.dev"]
beta.hi = high.income$coefficients[row.names(high.income$coefficients)=="avg.tmin.dev"]

# cluster robust standard errors
se.low = low.income$cse[names(low.income$cse)=="avg.tmin.dev"]
se.hi = high.income$cse[names(high.income$cse)=="avg.tmin.dev"]

# upper and lower confidence bounds
z.score.low = qnorm(1 - ((1 - 0.68)/2))
upper.low = beta.low + z.score.low*se.low
lower.low = beta.low - z.score.low*se.low
z.score.hi = qnorm(1 - ((1 - 0.68)/2))
upper.hi = beta.hi + z.score.hi*se.hi
lower.hi = beta.hi - z.score.hi*se.hi

d <- as.data.table(cbind(rbind(beta.low,beta.hi),rbind(upper.low,upper.hi),rbind(lower.low,lower.hi),rbind("Under $50,000","$50,000 and Over")))
setnames(d, c('beta', 'upper_bound', 'lower_bound', 'variables'))
d <- d[, beta := as.numeric(beta)*100] #rescale
d <- d[, upper_bound := as.numeric(upper_bound)*100]
d <- d[, lower_bound := as.numeric(lower_bound)*100]
d <- d[, variables := factor(variables, levels=c("Under $50,000", "$50,000 and Over"))]

inc <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        ggtitle("b") +
        geom_bar(stat = "identity", aes(width = 0.80), color='transparent') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("firebrick","gray70"),breaks=c("Under $50,000","$50,000 and Over")) +
        coord_cartesian(ylim = c(-1,10), expand=TRUE) +
        ylab("Marginal Effect of \u2103 Anomaly\non Monthly Nights of Insufficient Sleep\n(per 100 individuals)") +
        xlab("Income") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=20, angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=16, angle=0, vjust=.5, face="bold"),
              axis.text.x = element_text(size=15, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.y = element_line(),
              axis.line.x = element_line(),
              legend.position="none"
        )
 
# age
young <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
              | time + mmsa:season
              | 0 
              | mmsa + time
              , data=brfss[age < 65,], exactDOF=TRUE, psdef=FALSE)
summary(young)

old <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
            | time + mmsa:season
            | 0 
            | mmsa + time
            , data=brfss[age >= 65,], exactDOF=TRUE, psdef=FALSE)
summary(old)

# get coefficients of variables
beta.low = young$coefficients[row.names(young$coefficients)=="avg.tmin.dev"]
beta.hi = old$coefficients[row.names(old$coefficients)=="avg.tmin.dev"]

# cluster robust standard errors
se.low = young$cse[names(young$cse)=="avg.tmin.dev"]
se.hi = old$cse[names(old$cse)=="avg.tmin.dev"]

# upper and lower confidence bounds
z.score.low = qnorm(1 - ((1 - 0.68)/2))
upper.low = beta.low + z.score.low*se.low
lower.low = beta.low - z.score.low*se.low
z.score.hi = qnorm(1 - ((1 - 0.68)/2))
upper.hi = beta.hi + z.score.hi*se.hi
lower.hi = beta.hi - z.score.hi*se.hi

d <- as.data.table(cbind(rbind(beta.low,beta.hi),rbind(upper.low,upper.hi),rbind(lower.low,lower.hi),rbind("Under 65","65 and Over")))
setnames(d, c('beta', 'upper_bound', 'lower_bound', 'variables'))
d <- d[, beta := as.numeric(beta)*100] # rescale
d <- d[, upper_bound := as.numeric(upper_bound)*100]
d <- d[, lower_bound := as.numeric(lower_bound)*100]
d <- d[, variables := factor(variables, levels=c("Under 65","65 and Over"))]

age <- ggplot(data = d, aes(variables, beta, fill=variables)) +
        ggtitle("c") +
        geom_bar(stat = "identity", aes(width = 0.80), color='transparent') +
        geom_linerange(aes(ymax = upper_bound, ymin = lower_bound)) +
        scale_fill_manual(values=c("gray70","firebrick"),breaks=c("Under 65","65 and Over")) +
        coord_cartesian(ylim = c(-1,10), expand=TRUE) +
        ylab("Marginal Effect of \u2103 Anomaly\non Monthly Nights of Insufficient Sleep\n(per 100 individuals)") +
        xlab("Age") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=20,angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=16,angle=0, vjust=.5, face="bold") ,
              axis.text.x = element_text(size=15, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_line(),
              legend.position="none"
        )   

grid.arrange(season, inc, age, ncol = 3)
```

```{r echo=TRUE, eval=TRUE}
acute <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud  
               | mmsa + time 
               | 0 
               | mmsa + time
               , data=brfss[season=="summer" & income2 %in% 1:6 & age >= 65,], exactDOF=TRUE, psdef=FALSE)
summary(acute)
  
non.acute <- felm(qlrest2 ~ avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud  
               | mmsa + time 
               | 0 
               | mmsa + time
               , data=brfss[!(season=="summer" & income2 %in% 1:6 & age >= 65),], exactDOF=TRUE, psdef=FALSE)
summary(non.acute)
```

```{r figure3, echo=TRUE, fig.width=16, fig.height=6.5, dev='jpeg', dpi=1200, fig.cap="**Main text Figure 3**"}
# need to remove years from model names
bcsd <- bcsd[, model := gsub("2010", "", model)]
bcsd <- bcsd[, model := gsub("2050", "", model)]
bcsd <- bcsd[, model := gsub("2099", "", model)]

mmsa.2010 <- bcsd[year==2010, .(upper2010 = quantile(avg.tmin.dev, probs = 0.975, na.rm=T),
                    lower2010 = quantile(avg.tmin.dev, probs = 0.025, na.rm=T),
                    median2010 = quantile(avg.tmin.dev, probs = 0.5, na.rm=T),
                    mean2010 = mean(avg.tmin.dev, na.rm=T),
                    uppersem2010 = quantile(avg.tmin.dev, probs = 0.841, na.rm=T),
                    lowersem2010 = quantile(avg.tmin.dev, probs = 0.159, na.rm=T)), by=.(mmsa, model)]
mmsa.2050 <- bcsd[year==2050, .(upper2050 = quantile(avg.tmin.dev, probs = 0.975, na.rm=T),
                                lower2050 = quantile(avg.tmin.dev, probs = 0.025, na.rm=T),
                                median2050 = quantile(avg.tmin.dev, probs = 0.5, na.rm=T),
                                mean2050 = mean(avg.tmin.dev, na.rm=T),
                                uppersem2050 = quantile(avg.tmin.dev, probs = 0.841, na.rm=T),
                                lowersem2050 = quantile(avg.tmin.dev, probs = 0.159, na.rm=T)), by=.(mmsa, model)]
mmsa.2099 <- bcsd[year==2099, .(upper2099 = quantile(avg.tmin.dev, probs = 0.975, na.rm=T),
                                lower2099 = quantile(avg.tmin.dev, probs = 0.025, na.rm=T),
                                median2099 = quantile(avg.tmin.dev, probs = 0.5, na.rm=T),
                                mean2099 = mean(avg.tmin.dev, na.rm=T),
                                uppersem2099 = quantile(avg.tmin.dev, probs = 0.841, na.rm=T),
                                lowersem2099 = quantile(avg.tmin.dev, probs = 0.159, na.rm=T)), by=.(mmsa, model)]
setkey(mmsa.2010, mmsa, model)
setkey(mmsa.2050, mmsa, model)
setkey(mmsa.2099, mmsa, model)

mmsa.bcsd <- merge(mmsa.2010, mmsa.2050)
mmsa.bcsd <- merge(mmsa.bcsd, mmsa.2099)
rm(mmsa.2010,mmsa.2050,mmsa.2099)

# get tmin mean differences due to climate change between 2010 and 2050 and 2050 and 2099
mmsa.bcsd <- mmsa.bcsd[, ':=' (mean.dif.2050.2010 = mean2050-mean2010)]
mmsa.bcsd <- mmsa.bcsd[, ':=' (mean.dif.2099.2010 = (mean2099-mean2050) + mean.dif.2050.2010)] # incremental difference from 2010 to 2099 for plotting.
mmsa.bcsd <- mmsa.bcsd[, baseline := 0]

# line plot data
lineplot <- mmsa.bcsd[, .(mmsa = mmsa, model=model, baseline = baseline, diff2050 = mean.dif.2050.2010, diff2099 = mean.dif.2099.2010)]
lineplot <- melt(lineplot, id.vars=c("mmsa","model"))
setnames(lineplot, c("variable","value"), c("year", "clim.tmin"))
lineplot <- lineplot[year=="baseline", year := "2010"]
lineplot <- lineplot[year=="diff2050", year := "2050"]
lineplot <- lineplot[year=="diff2099", year := "2099"]
lineplot <- lineplot[, year := as.numeric(as.character(year))]
lineplot <- lineplot[, id := .GRP, by=.(mmsa, model)]
setkey(lineplot, id, year)

# create forecasts for 2050 and 2099 of monthly nights of sleep lost per 100 individuals
model <- felm(qlrest2 ~  avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
              | time + mmsa:season
              | 0 
              | mmsa + time, data=brfss, exactDOF=TRUE, psdef=FALSE)
tmin.coef <- model$coefficients[1]
lineplot <- lineplot[, forecast := clim.tmin*tmin.coef*100]
lineplot2050 <- lineplot[year==2010 | year==2050, ]
lineplot2099 <- lineplot[year==2050 | year==2099, ]

# data for median forecast line
med.line2050 <- lineplot2050[, .(forecast=median(forecast, na.rm=T)), by=.(year)]
med.line2099 <- lineplot2099[, .(forecast=median(forecast, na.rm=T)), by=.(year)]

# data for vertical lines
vline2050 <- lineplot[year==2050, .(min=min(forecast), max=max(forecast), year=2050)]
vline2099 <- lineplot[year==2099, .(min=min(forecast), max=max(forecast), year=2099)]

lineforecast <- ggplot(data=NULL, aes(x=year, y=forecast, group=id)) +
    geom_line(data=lineplot2050, color="goldenrod1", alpha=0.02) +
    geom_line(data=med.line2050, aes(group=NULL), linetype="dashed") +
    geom_line(data=lineplot2099, color="firebrick", alpha=0.02) +
    geom_line(data=med.line2099, aes(group=NULL), linetype="dashed") +
    geom_linerange(data=vline2050, aes(x=year, ymin=min, ymax=max, y=NULL, group=NULL)) +
    geom_linerange(data=vline2099, aes(x=year, ymin=min, ymax=max, y=NULL, group=NULL)) +
    scale_x_continuous(breaks=c(2010,2050,2099)) +
    coord_cartesian(ylim = c(-0.0025,23), expand=TRUE) +
    ggtitle('b') +
    ylab("Predicted Monthly Nights of Insufficient Sleep\n(per 100 Individuals)") +
    xlab("Forecast Year\n(Baseline 2010)") +
    theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.025),
          axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=17,angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
          axis.text.x = element_text(size=18, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank()
    )

# bar plot (mean and SEM) data
bcsd <- bcsd[, forecast := avg.tmin.dev*tmin.coef*100]
barplot <- bcsd[, .(mean = mean(forecast, na.rm=T),
                    uppersem = quantile(forecast, probs = 0.841, na.rm=T),
                    lowersem = quantile(forecast, probs = 0.159, na.rm=T)), by=.(year)]
barplot <- barplot[year!=2010,]
barplot <- barplot[, year := factor(year, levels=c("2050", "2099"))]

barforecast <- ggplot(data = barplot, aes(x=year, y=mean, fill=year)) +
        ggtitle("c") +
        geom_bar(stat = "identity", aes(width = 0.8), color='transparent') +
        geom_linerange(aes(ymax = uppersem, ymin = lowersem)) +
        scale_fill_manual(values=c("goldenrod1","firebrick"),breaks=c("2050","2099")) +
        coord_cartesian(ylim = c(-0.0025,23), expand=TRUE) +
        ylab("Mean Predicted Monthly Nights of Insufficient Sleep\n(per 100 Individuals)") +
        xlab("Forecast Year\n(Baseline 2010)") +
        geom_hline(yintercept = 0) +
        theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
              plot.title = element_text(hjust=0.025),
              axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
              axis.title.y = element_text(size=17, angle=90, vjust=1, face="bold"),
              axis.text.y = element_text(size=18, angle=0, vjust=.5, face="bold"),
              axis.text.x = element_text(size=18, face="bold", color="gray30"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = element_line(),
              legend.position="none"
        )

# histograms of nightly temperature deviations from baseline
histforecast <- ggplot() +
                  geom_histogram(data=bcsd[year==2010,], aes(x=avg.tmin.dev, y = ..density.., fill="gray40"), color="gray30", alpha=.75) +
                  geom_histogram(data=bcsd[year==2050,], aes(x=avg.tmin.dev, y = ..density.., fill="goldenrod1"), color="goldenrod4", alpha=.75) +
                  geom_histogram(data=bcsd[year==2099,], aes(x=avg.tmin.dev, y = ..density.., fill="firebrick"), color="firebrick4", alpha=.75) +
                  ylab("Density") +
                  xlab("Nightly Temperature\n(Anomaly in \u2103)") +
                  coord_cartesian(ylim=c(-0.0025, 0.31), expand=TRUE) +
                  scale_x_continuous(breaks=c(-10,-5,0,5,10,15)) +
                  ggtitle("a") +
                  scale_fill_manual(limits=c("gray40", "goldenrod1", "firebrick"), values=c("gray40", "goldenrod1", "firebrick"),labels=c("2010", "2050", "2099")) +
                  guides(fill = guide_legend(override.aes = list(colour = NULL))) +
                  theme(title = element_text(size=20, angle=0, vjust=.5, face="bold"),
                        plot.title = element_text(hjust=0.025),
                        axis.title.x = element_text(size=20, angle=0, vjust=-0, face="bold"),
                        axis.title.y = element_text(size=17,angle=90, vjust=1, face="bold"),
                        axis.text.y = element_text(size=18,angle=0, vjust=.5, face="bold") ,
                        axis.text.x = element_text(size=18, face="bold"),
                        panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(),
                        axis.line.x = element_line(),
                        axis.line.y = element_line(),
                        legend.title=element_blank(),
                        legend.position = c(.15, 0.2),
                        legend.text=element_text(size=16, face="bold"),
                        legend.key = element_rect(colour = "black")
                  )

grid.arrange(histforecast, lineforecast, barforecast, ncol = 3)
```

```{r figure4, echo=TRUE, fig.width=8, fig.height=4, dev='jpeg', dpi=1200, fig.cap="**Main text Figure 4**"}
# get state shapefiles
setwd("./shp/")
us <- readOGR(dsn=getwd(), layer="cb_2014_us_state_5m", verbose=FALSE) #read in each spatial polygon data frame
us <- us[!(us$NAME=="Puerto Rico" | us$NAME=="United States Virgin Islands" | us$NAME=="American Samoa" | us$NAME=="Guam" | us$NAME=="Commonwealth of the Northern Mariana Islands" | us$NAME=="Alaska" | us$NAME=="Hawaii"),] #continental us
projection(us) <- projection(rast2050)

model <- felm(qlrest2 ~  avg.tmin.dev + I(avg.tmax - avg.tmin) + sum.prcp.dev + avg.humid + avg.cloud 
              | time + mmsa:season
              | 0 
              | mmsa + time, data=brfss, exactDOF=TRUE, psdef=FALSE)
summary(model)

tmin.coef <- model$coefficients[1]
forecast2050 <- rast2050*tmin.coef*100
forecast2099 <- rast2099*tmin.coef*100

# trim forecast rasters to remove outer NA values
forecast2050 <- trim(forecast2050, padding=0, values=NA)
forecast2099 <- trim(forecast2099, padding=0, values=NA)

# stack these together
forecast <- stack(forecast2050, forecast2099)
id <- c('2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2050','2099') # set names

# plot raster forecasts
theme = rasterTheme(region=c('gray40', 'goldenrod1', 'firebrick4')) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               layout=c(2,1), # horizontal layout
               names.attr=c("2050","2099"), # inset titles
               colorkey=list(space="bottom"), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
p <- update(p, par.settings = list(axis.line = list(col = 0)), sub="Predicted Monthly Nights of Insufficient Sleep\n(per 100 individuals)") # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(us, lwd=0.8, col="black")) # add us shapefiles

p
```
